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The  Marine  Gas  Turbine  control  systems  in  present  use  in  the  U.  S. 
Navy  are  of  such  significant  technological  age  that  new  design 
techniques  could  lead  to  more  optimum  performance  and  increased 
plant  efficiency.  To  this  end.  a  new  real  time  Marine  Gas  Turbine 
simulation  method  is  needed  for  advanced  controller  design  and 
implementation. 

A  modeling  method  is  shown  which  utilizes  real  time  sequential 
linearizations  to  approximate  the  true  nonlinear  response  of  the  NPS 
Boeing  502-6A  test  facility.  A  validation  of  this  simulation  approach  is 
presented.  The  method  has  immediate  application  to  advanced 
controller  design,  especially  to  the  design  of  modem  regulators 
(Linear  Quadratic),  model  reference  controllers,  and  real  time 
diagnostics. 
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E  -  Fuel  Combustion  Energy  Realized  at  HP  Turbine 
JD  -  Dynamometer  Torque 
JG  -  Gas  Generator  Torque 
Ma  -  Mass  Flowrate  of  Air 
Maf  -  Mass  Flowrate  of  Air/Fuel  Mixture 
Mf  -  Mass  Flowrate  of  Fuel 
Ng  -  Gas  Generator  Speed 
Ns  -  Power  Turbine/Dynamometer  Speed 
Qc  -  Compressor  Torque 
Qd  -  Dynamometer  Torque 
Qfpt  -  Free  Power  Turbine  Torque 
Qhpt  -  High  Pressure  Turbine  Torque 
P2  -  Compressor  Discharge  Pressure 
P4  -  High  Pressure  Turbine  Discharge  Pressure 
T2  -  Compressor  Discharge  Temperature 
T4  _  High  Pressure  Turbine  Discharge  Temperature 

Note:  1.  Lower  case  variables  represent  normalized 
variables. 

2.  _ B  denotes  normalizing  values.  ( i.e.  JGB,  NGB, 

QCB,  etc.) 
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Conventional  control  theory  is  limited  to  time  invariant  single -input- 
single-output  systems  and  is  usually  based  on  a  frequency  domain 
design  approach.  Modem  control  theory,  on  the  other  hand,  can  be 
well  applied  to  multiple-input-multiple-output,  non-linear  gas  turbine 
systems  which  are  now  in  use  and  it  relies  on  a  much  more  applicable 

time  doirain  approach.  1  The  feasibility  of  modern  control  theory 

such  as  Linear  Quadratic  Regulator  (LQR)  theory,  in  the  application  of 

marine  gas  turbines  for  propulsion  has  been  demonstrated  by  past 

work  at  the  Naval  Postgraduate  School  in  Monterey,  California. 2  4  For 

this  application  to  be  successful  a  simple,  accurate,  and  robust  real 

time  marine  gas  turbine  simulation  must  be  developed. 

This  thesis  introduces  a  modeling  approach  to  be  applied  in 
developing  a  real  time  simulation  of  engine  response  to  meet  the 
requirements  of  simplicity,  accuracy,  and  robustness.  It  is  simple  in 
that  it  is  only  first  order  dependent,  accurate  in  that  it  provides 
acceptable  (+/-  10%)  response  for  use  in  application  of  advanced 
engine  control  design,  and  robust  in  that  it  wo^ks  well  over  the  entire 
operational  range  of  the  gas  turbine  system. 


n.  PAST  WORK  IN  GAS  TURBINE  MODELING 

Past  work  in  Gas  Turbine  Modeling  is  discussed  in  terms  of  two 
areas:  first,  a  brief  discussion  on  gas  turbine  modeling  in  general  is 
presented.  This  is  followed  by  a  discussion  on  the  previous  work 
completed  at  the  Naval  Postgraduate  School  in  Monterey,  California. 

Since  most  marine  gas  turbines  have  been  derived  from  aviation 
sources,  most  of  the  early  gas  turbine  modeling  work  was 
accomplished  in  the  aviation  arena.  In  the  1970’s  the  Naval  Gas 
Turbine  Ship  Propulsion  Dynamics  and  Control  Systems  Research  and 
Development  Program  was  started  to  address  the  issues  of  designing 

Marine  Gas  Turbines  for  future  use  in  the  Navy. 5  This  included  the 

dynamic  modeling  of  the  response  of  j.  marine  gas  turbine.  Prior  to 

this  program,  most  of  the  work  had  been  steady  state  modeling. 

However,  the  Navy  realized  that  dynamic  modeling  was  needed  for 

advanced  controller  design.  The  program  was  successful  in  developing 

a  gas  turbine  machinery  dynamics  and  control  system  data  base.  From 

this  work  a  computer-based  propulsion  control  testbed  was  developed 

for  use  in  gas  turbine  controller  development.  Cut  and  try  linear 

control  design  then  ensued,  which  resulted  in  a  gain  adjusted 
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controller  design  in  use  in  today’s  fleet.  This  approach  resulted  in  a 
linear  Proportional-Integral  theory  being  applied  to  a  nonlinear 
system.  The  results  were  deemed  adequate  for  the  Navy's  uses. 

A  new  approach  to  this  modeling  is  to  develop  a  method  .o 
accurately  model  the  propulsion  plant  nonlinear  responses  with  a  real 
time  simulator.  This  would  help  eliminate  or  alleviate  problems 
associated  with  the  varying  propeller  loads  felt  in  the  gas  turbine  by 
predicting  the  responses  before  they  actually  occur,  thus  allowing 
more  acurate  control.  More  accurate  control  would  mean  better 
engine  response  to  these  loads  and  would,  in  turn,  mean  less  wear, 
thus  leading  to  longer  engine  life  as  well  as  better  engine  fuel 
efficiency.  This  ir.  the  direction  taken  at  the  Naval  Postgraduate 
School,  (NPS),  in  the  gas  turbine  modeling  problem. 

The  previous  work  started  with  Phillip  Johnson  in  1985  in  which 
he  conducted  the  development  and  implementation  of  a  load  control 
system  (  e.g.  a  water  brake  dynamometer)  to  emulate  the  shaft  and 
propeller  of  a  marine  gas  turbine.®  The  nexi  step  was  taken  by  James 
Roger  also  in  1985.  He  used  a  simple  multiport  diagram  and  the 
Continous  System  Modeling  Program  (CSMP  III)  to  incorporate  an 
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improved  model  for  the  dynamometer  unloading  valve  into  the  existing 

system  dynamic  model. ?  These  two  works  formed  the  base  for  the 
next  step  performed  by  Vincent  Herda  in  1986.  Herda  developed  a 
better  multiport  model  and  used  it  to  develop  a  linear  state  space 

model  and  a  nonlinear  dynamic  propulsion  model.2  Robert  Miller 
followed  in  the  same  year  using  Herda's  models  and  examined  the 
feasibility  of  using  Linear  Quadratic  Regulator  (LQR)  control  design 
theory  for  a  small  gas  turbine,  specifically  the  Boeing  502-6A  test  bed 

installation  at  NPS.3  Finally  Vincent  Stammetti  in  1988  developed  the 
first  real  time  dynamic  simulator  for  the  gas  turbine  and  compared  an 
LQR  controller  to  a  classical  Proportional  Integral  (PI)  regulator  based 

on  his  results  with  the  simulator.4  To  summarize.  Herda  validated 
cause  and  effect  through  the  use  of  his  new  multiport  model,  shown  in 
figure  1  ,  Miller  validated  the  steady-state  model  method  and 
Stammetti  followed  with  the  first  linear  model  structure.  It  is 
important  to  note  that  these  models  were  only  accurate  for  one  run 
of  simulation  and  were  not  applicable  over  the  whole  range  ol 
operation  of  the  gas  turbine.  Their  main  point  however,  was  to 
validate  and  show  the  feasibility  of  developing  a  better  model  for  use  in 
modern  controller  design  to  update  those  in  current  use  today.  To 
this  end  they  were  successful.  This  paper  will  continue  this  work  and 
attempt  to  develop  a  simple,  robust,  and  accurate  real  time  marine  gas 


4 


:erda’s  Multiport  Diagram. 


m.  MODEL  STRUCTURE  AND  COARSE  TUNING 

A.  Modeling 

The  gas  turbine  plant  to  be  modeled  was  the  175  horsepower 
Boeing  502 -6A  which  was  assembled  in  a  test  bed  configuration  at  the 
Naval  Postgraduate  School  in  Monterey,  California.  It  was  coupled  to  a 
Clayton  17-300  water  brake  dynamometer  which  was  used  to  simulate 
the  shaft  and  propeller  loading  on  a  marine  gas  turbine.  The  gas 
turbine /water  brake  assembly  is  shown  in  figure  2.  The  gas  turbine 
itself  was  composed  of  two  aerodynamically  coupled  sections,  a  gas 
generating  section  and  a  power  output  section.  The  gas  generator  had 
three  major  components.  These  were  the  compressor,  the  burner, 
and  the  high  pressure  turbine  shown  in  figure  3.  Connected  by 
gearing  off  the  gas  generator  shaft  was  the  accessory  section 
consisting  of  the  fuel  pump,  lube  oil  pump,  governor,  tachometer 
generator,  and  the  starter.  The  compressor  was  a  single  stage 
centrifugal  compressor  and  was  coupled  to  the  high  pressure  turbine 
aerodymanically  via  two  through-flow  type,  cross-connected 
combustion  chambers.  The  high  pressure  turbine  was  a  single  stage 
axial  flow  high  pressure  turbine.  It  was  aerodynamically  connected  to 
the  low  pressure  or  free  power  turbine  which  was  also  an  axial  flow 
turbine.  The  free  power  turbine  output  shaft  was  mechanically 
coupled  to  the  water  brake  dynamometer  which  absorbs  the  energy  of 
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175  Hp  Boeing  502-6A  gas  turbine 

Single  stage  centrifugal  compressor 
Two  single  stage  axial  turbines 


Clayton  17-300  water  brake  dynamometer 


Assorted  sensors 
Pressure 
Temperature 
Output  torque 
Rotational  speeds 
Fuel  flow  rate 
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Figure  3.  Gas  Turbine  Components. 


the  free  power  turbine  in  simulating  the  action  of  a  drive  shaft  and 
propeller. 

The  modeling  began  with  Herda’s  cause  and  effect  multiport 
diagram  shown  in  figure  1.  The  importance  of  this  diagram  to  the 
method  can  hardly  be  overstated.  Each  of  the  components,  the 
compressor,  high  pressure  and  free  power  turbines,  and  the  gas 
generator  and  dynamometer  shafts,  were  modeled  by  describing  their 
outputs  in  terms  of  their  inputs  in  equation  form.  For  example,  the 
compressor  inputs  are  gas  generator  speed  (Ng)  and  the  high  pressure 
turbine  inlet  pressure  (P2)  and  the  outputs  are  mass  flow  rate  of  air 
(Ma),  high  pressure  turbine  inlet  temperature  (T2),  and  compressor 
torque  (Qc).  The  compressor  block  diagram  is  shown  in  figure  4. 

Ma 
T2 

P2  Figure  4.  Compressor  Block  Diagram. 

It  is  important  to  note  that  these  input/output  quantities  as  well  as 
those  for  the  other  components  must  be  measurable  or  calculable. 
They  were  recorded  and  or  calculated  over  the  operating  range  of  the 
gas  turbine  (shown  in  figure  5),  and  the  results  are  tabulated  in  table 
1.  For  clarity,  the  operating  region  of  the  turbine  is  shown  in  figure  5. 
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Figure  5.  Gas  Turbine  Operating  Range. 
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0.88538  1.67267  0.18182  0.81386  0.66643  0.93828  0.68421  0.77535  0  85169  0.90542  0.84148 

1.00038  0.33733  1.11364  0.95767  0.87143  0.92678  0.83659  0.88416  1.02542  0.»103S  1.01243 

1.00346  0.67067  0.90909  0.96084  0.89643  0.96444  0.90526  0.90473  1.00047  0.89222  0.99619 

1.00308  1.00133  0.77273  0.96758  0.88786  0.95920  0.94737  0.92735  0.99576  0.92797  0.98411 

1.00346  1  33333  0.59091  0.96465  0.89643  0.96339  1.00000  0.92735  0.99576  0.95382  0.98411 

1.00346  1.67333  0.45455  0.96977  0.89643  0.96862  0.97368  0.92735  1.00000  0.94584  0.°?825 

1.11615  0.33933  1.4545S  1.13572  1.14643  0.«5502  1.10526  1.05037  1.16102  1.07252  1.14700 

1.11769  0.67333  1.27273  1.15037  1.14643  1.00418  1.15789  1.07298  1.15"8S  1.00119  1.12666 

1.11577  1.00267  1.04545  1.14349  1.16786  1.01151  1.26316  1.10588  1.13136  0.48389  1.11886 

1.11731  1.33333  0.88636  1.13656  1.17500  1.01464  1.31579  1.11404  1.12712  1.01459  1.11484 

1.11769  1.67667  0.72727  1.12744  1.17143  1.00000  1.31579  1.12033  1.12712  1.07170  1.11493 

1.23192  0.33933  1.90909  1.30535  1.46428  1.02301  1.42105  1.28849  1.29661  1.24762  1.28258 

1.23154  0.67333  1.65909  1.31642  1.48571  1.04812  1.50000  1.32131  1.283»0  1.18356  1.27065 

1.23154  1.00667  1.40909  1.30721  1.49286  1 . 086B2  1.57895  1.34784  1.27542  1.12776  1.26275 

1.23115  1.34467  1.25000  1.31860  1.49643  1.07113  1.73684  1.36816  1.28814  1.14542  1.27544 

1.23077  1.67333  1.09091  1.32419  1.48928  1.0«205  1.73684  1.37420  1.26695  1.12057  1.25486 

(B) 


(A)  Raw  data 
(B)  Normalized  data 

Table  1 .  Steady  State  Data  Taken  Over  Operating  Range 
(Ng  20k  -32k  rpm,  Ns  500  -  2500  rpm) 
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At  this  point  it  should  be  noted  that  the  data  of  table  1A  was 
normalized  according  to  each  of  the  variable’s  midpoint  for  ease  of 
understanding  each  variable's  relative  impact  with  respect  to  the 
others.  Once  normalized,  the  variables  were  are  denoted  using  lower 
case  letters.  The  resulting  equation  forms  are: 


Ma  =  fi  (Ng,  P2), 

(1) 

T2  =  f2  (Ng,  P2), 

and 

(2) 

Qc  =  f3  (Ng,  P2). 

(3) 

The  normalization  values  are  shown  in  table  2. 

Table  2.  Normalization  Points. 

Ng  =  26000  Ns  =  1500 

Qd/Qfpt  =  220  T2  =  215 

P2  =  14  T4  =  956 

P4  =1.9  Mf  =  122.5 

Ma  =  2.36  Qc/Qhpt  =  80 

The  next  step  was  to  distinguish  between  steady  and  dynamic 
components.  Once  this  is  complete,  the  component  modeling  can 
proceed.  The  rest  of  the  component  equation  forms  are  shown  in 
table  3. 


Table  3.  Component  Equations. 


Gas  Generator  Shaft: 

Ng  = 

=  U  (Qhpt,  Qc) 

(4) 

Dynamometer  Shaft: 

Ns  = 

fs  (Qfpt.  Qd) 

(5) 

Free  Power  Turbine: 

P4  = 

fe  (Ns,  T4,  Maf) 

(6) 

Qfpt  = 

{7  (Ns,  T4,  Maf) 

(7) 

High  Pressure  Turbine:  E  = 

fs  (Mf) 

(8) 

Maf  = 

f9  (Ma,  Mf) 

(9) 

T4  = 

fio  (Ma.  T2,  P4, 

Ng.  E) 

(10) 

P2  = 

fi  1  (Ma,  T2,  P4, 

Ng,  E) 

(ID 

Qhpt  = 

to 

S 

ft) 

s 

Ng.  E) 

(12) 

B.  State  Selections 

State  variables  are  the  smallest  set  of  variables  which  are  required 
to  ascertain  the  state  or  condition  of  a  dynamic  system.1  The 
knowledge  of  these  variables  at  an  initial  point  in  time  (to),  when 

combined  with  the  input  as  time  proceeds  (t  >  to),  will  completly 
describe  the  systems  behavior.  The  state  of  a  system  at  any  time  (t)  is 
indepentent  of  any  prior  state  and  input  which  precede  the  initial 

time  point  (to) 1 .  Here,  the  turbine  and  load  is  the  dynamic  system  and 
the  states,  when  properly  applied,  will  accurately  describe  the 


dynamic  behavior  of  the  system.  Significant  prior  research  has 
determined  that  only  three  states  are  necessary  to  sufficiently 

describe  the  operating  condition  of  the  NPS  gas  turbine.®  In  reality 
however,  only  two  states  are  necessary  (ng,  and  ns)  since  the  third 
state  (e)  represents  a  short  first  order  lag  of  the  input  fuel  (mf).  This 
time  lag  is  due  to  the  burner  action  delay  between  the  input  of 
additional  fuel  and  the  thermodynamic  energy  output  felt  in  the 
turbine. 

It  is  our  goal  to  relate  all  of  the  dependent  variables  (p2. 
p4,t2.mf,...)  to  the  three  state  variables  ng,  ns.  and  e  to  allow  us  to  use 
the  state  space  equation  x  =  Ax  +  Bu  (ref.  9)  to  describe  the  transient 
behavior  of  our  system.  Here  x  is  the  state  vector  and  is  composed  of 
state  variables  ng,  ns,  and  e.  u  is  the  input  vector  and  is  composed  of 
mf  and  qd.  Refering  back  to  figure  1,  note  that  we  have  redrawn 
Herda's  system  boundary  at  qd  in  order  to  have  a  more  measurable 
input.  Again,  from  prior  research,  it  was  known  that  each  of  the  states 
is  associated  with  a  significant  time  lag  in  the  system  response.  In 
this  way,  we  knew  that  the  fuel  input  component,  the  rotor  shaft,  and 
the  output  shaft  must  be  regarded  as  dynamic  components.  We  also 
knew  that  the  remaining  components  could  be  regarded  as 
instantaneous  (memoryless),  or  steady  components.  The  steady  and 
dynamic  components  are  shown  in  table  4. 


Table  4.  Steady  and  Dynamic  Components. 

Steady  Components: 

compressor  ma 

high  pressure  turbine  p4 
power  turbine  t4 

Dymamic  Components: 

rotor  shaft  ng 

load  shaft  ns 

burner  lag  e 

C.  Steady  Components 

While  the  steady  components  of  the  model  were  given  above  as  non 
linear  functions,  in  this  section,  we  reduce  these  to  linear 
expressions,  where  the  coefficients  were  regressed  and  varried 
according  to  the  engine  operating  point. 

The  steady  component  equation  form  selection  was  accomplished 
through  the  use  of  an  equation  fitting  program  (Minitab),  and  plotting 
of  each  of  the  input  variables  versus  an  output  variable  minus  the 
residual.  If  the  residual  plot  displayed  a  parabolic  shape,  it  was  an 
indication  that  the  input  variable  function  in  question  was  of  a  higher 
order.  In  the  case  of  the  equations  with  only  two  input  variables  a 
different  method  was  used  to  fit  the  data.  That  is.  one  at  a  time,  each 


t2  qc 
qfpt 

p2  qhpt  maf 


of  the  input  variables  were  held  constant  and  the  other  was  plotted 
versus  the  output  variable.  Both  methods  are  developed  further  in 
Appendix  C.  This  resulted  in  accurate  data  fit  equations. 

To  simplify  the  equation  reduction  methodology,  the  steady 


component  equations  were  kept  in  the  form: 

x  =  cl  y  +  c2  z  (13) 

For  example,  the  compressor  equations  appeared  as  follows: 

t2  =  cl  ng2  +  c2  ng  +  c3  p2  +  c4  (14) 

ma  =  c5  ng2  +  c6  ng  +  c7  p2  +  c8  and  (15) 

qc  =  c9  ng2  +  clO  ng  +  cl  1  p2  +  cl2  (16) 

where  the  coefficients  cl  thru  cl2  were  the  values  from  the  equation 


fits.  It  is  important  to  note  the  strong  dependency  of  the  state  ng. 
This  will  appear  again  in  the  decoupling  equations  developed  at  the 
end  of  this  chapter.  The  next  phase  of  development  was  the 
formulation  of  the  dynamic  component  models. 

D.  Dynamic  Components 

The  dynamic  components  are  derived  from  the  gas  generator  shaft, 
the  dynamometer  shaft,  and  the  lag  relationship  between  mf  and  e. 
The  dynamic  components  were  modeled  with  differential  equations 
which,  when  integrated  over  time  will  determine  the  state  of  the 


system.  The  dynamic  components  took  the  form: 


ng  =  qhpt-qc 

(17) 

ns  =  qfpt-qd 

(18) 

and 

e  =  (mf-e)/T. 

(19) 

These,  along  with  the  previous  equations  of  the  steady  components, 
were  next  linearized  to  form  a  dynamic  equation  set  for  the  simulation 
model. 

E.  Linearized  Dynamic  Equation  Set 

The  linearized  dynamic  equation  set  was  developed  from  the 
perturbational  or  derivative  values  of  the  steady  state  equations 
(symbolically,  the  perturbational  variables  will  be  proceeded  by  a  "d"). 
The  compressor  equations  are  used  below  to  demonstrate  the 
formulation  of  the  dynamic  equations.  Recall  that  the  t2  steady  state 
equation  was: 

t2  =  cl  ng2  +  c2  ng  +  c3  p2  +  c4  (14) 

So.  the  dynamic  equation  took  the  form: 

dt2  =  (  dt2./')ng  j  dng  +  [  8t2/9p2  ]  dp2  (20) 

Here  the  partial  differentials  will  be  denoted  with  upper  case  letters 
so  that  the  resulting  form  of  equation  20  is: 

dt2  =  Al  dng  +  A2  dp2. 
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(21) 


The  rest  of  the  steady  and  dynamic  equations  previously  developed 
were  also  converted  to  dynamic  form  in  this  manner.  The  results 
appear  below  in  table  5. 

Table  5.  Dynamic  Equation  Set. 


Compressor  dt2  =  Al  dng  +  A2  dp2  . (22) 

dma  =  B  l  dng  +  B2  dp2  . (23) 

dqc  =  Cl  dng  +  C2  dp2  .  (24) 

Gas  Generator  Shaft 

JGB  dng  =  dqhpt-dqc . (25) 

Dynamometer  Shaft 

JDB  dns  =  dqfpt-dqd  . (26) 

Free  Power  Turbine 

dp4  =  D1  dns+D2  dt4+D3  dmaf .  (27) 

dqfpt  *Fl  dns+F2  dt4+F3  dmaf  .  (28) 

High  Pressure  Turbine 

de  =  (dmf-de)/T  . (29) 

from  continuity—  dmaf  =  Kl  dmf  +  K2  dma .  (30) 

dt4  =  Gl  dma  +  G2  dt2  +  G3  dp4  +  G4  dng  +  G5  de  .  (31) 

dp2  =  Hi  dma  +  H2  dt2  +  H3  dp4  +  H4  dng  +  H5  de  .  (32) 

dqhpt  =  U  dma  +  12  dt2  +  13  dp4  +  14  dng  +  15  de  .  (33) 


These  resulting  equations  from  table  5  will  be  used  to  develop  the 
state  space  relationships  through  a  reduction  methodology  to  get 
three  equations  (  one  each  for  ng,  ns.  and  e)  in  terms  of  ng,  ns.  e.  mf, 


and  qd. 


F.  Decoupling  the  Dynamic  Equation  Set 

In  the  past,  during  the  equation  set  reduction,  a  number  of  the  table 
3.  equations  were  used  more  than  once.  This  resulted  in  the 
development  of  singularity  points  in  the  simulations.  To  avoid  this, 
another  method  was  developed  in  the  present  work.  This  was 
achieved  by  inspecting  the  system  block  diagram  (fig.  1 1  and 
deciding  that  if  two  equations,  one  for  p2  and  the  other  for  p4, 
could  be  expressed  in  terms  of  the  states  ng,  and  ns  only,  then  it 
would  be  possible  to  reduce  the  equation  set  directly  (without 
resubstitution)  thus  eliminating  the  singularities.  The  new  p2  and  p4 
equations  weie  developed  with  the  assistance  of  the  Minitab 
regression  routine  and  steady  state  data  for  ng  and  ns.  This  resulted 
in  the  conclusion  that  both  variables  had  a  second  power  effect  with 
ng  and  a  first  power  effect  with  ns.  With  these  relationships  in  mind, 
the  computer  routine  was  used  to  fit  the  data  in  a  two  step  process. 

First  ,  equations  for  p2  and  p4  were  generated  in  terms  of  ngA2,  ng, 
and  ns.  Next,  the  data  generated  from  these  equations  for  p2  and  p4 
was  plotted  versus  ng  at  five  different  constant  ns  values.  These  plots 
were  compared  to  the  actual  data  to  check  for  accuracy  of  the 
equations.  The  results  are  shown  on  the  two  graphs  in  figure  6  and 
show  the  robustness  of  these  equations  over  the  entire  range  of  the 
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P2  VS  NG  (NS=500-2500  RPM) 


o 


20000  23000  26000  20000  32000 

NG  (RPM) 


P4  VS  NG  (NS=500— 2500  RPM) 


*0000  23000  20000  29000  32000 


NG  (R PM) 

Figure  6.  P2  and  P4  graphical  Representations. 
Note:  Solid  lines  are  equations,  symbols  are  data. 
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data.  The  new  dynamic  equations  for  p2  and  p4  are  now  in  the  form  of: 


p2  =  cl  ngA2  +  c2  ng  +  c3  ns  +  c4  and  (34) 

p4  =  c5  ngA2  +  c6  ng  +  c7  ns  +  c8.  (35) 

These  equations  were  then  linearized  into  the  form: 

dp2  =  Hi  dng  +  H2  dns  for  eqn.  11  and  (36) 
dp4  =  Dl  dng  +  D2  dns  for  eqn.  6.  (37) 

These  twelve  equations  were  then  reduced  with  the  reduction 
methodology  as  described  in  detail  in  Appendices  A  and  B. 

The  results  of  the  reduction  transformed  the  linearized  equations  to 
state  space  equations,  which  were  used  in  the  simulation  model.  The 
resulting  state  space  equations  from  Appendix  B  are: 

dng  =  Z1  dng  +  Z2  dns  +  Z3  de 
JGB  JGB  JGB 

dns  =  Z4  dng  +  dns  +  Zfi  de  +  Z7  dmf  -  1  dqd 
JDB  JDB  JDB  JDB  JDB 

de  =  -  1  de  +  1  dmf 

T  T 

where  Z1  through  Z7  were  derived  from  combinations  of  the  A1 
through  K2  coefficients  of  the  dynamic  equations.  The  standard  state 

space  equation  is  x  =  A(t)  x  +E.(t)  u.9  Here  &  is  a  matrix  made  up  of 
dng,  dns,  and  de,  x  is  a  state  vector  composed  of  the  state  variables: 
dng,  dns,  and  de,  and  u  is  a  vector  composed  of  inputs  dmf  and  dqd. 
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The  A  and  B  matrices  are: 

Al  1  =  Zi  A12  =  Z2  Al3  =  Z2. 

JGB  JGB  JGB 

A21=Z4  A22  =  ZZl  A23  =  Z3. 

JDB  JDB  JDB 

A31  =  O  A32  =  O  A3  3  =  zl_ 

T 

and  Bll  =  O  Bl2  =  O 

B2 1  =  Z7  B22  =  -1 

JDB  JDB 

B31  =  _1_  B32  =  O 

T 

and  the  state  space  equations  now  appear  as: 

dng  =  All  dng  +  Al2  dns  +  A13  de  +  Bl  1  dmf  +  Bl2  dmf 
dns  =  A2 1  dng  +  A22  dns  +  A23  de  +  B2 1  dmf  +  B22  dmf 

de  =  A31  dng  +  A32  dns  +  A3 3  de  +  B31  dmf  +  B32  dmf. 

A  coarse  tuning  was  next  performed  to  determine  the  best 
numerical  values  of  the  coefficients.  This  was  necessary  since  the  first 
set  of  regressed  coefficients  did  not  provide  accurate  simulation.  This 
was  attributed  to  the  coarse  instrumentation  that  is  present  on  the 
NPS  turbine  and  the  manner  of  fitting  in  the  regression  routine. 
Consequently,  attention  was  focused  on  the  following  coefficients  as  in 
need  of  tuning:  Al,  A2,  Bl,  B2,  Cl,  C2,  FI.  and  II  -  15  .  Those  for  Dl, 
D2,  HI,  H2,  Kl,  and  K2  were  found  previously  using  the  data 
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regression  program.  Initially  chosen  for  a  starting  point  was  that  for 
B2  (9ma/9p2).  This  was  an  obvious  choice  for  much  is  known  about 
the  relationship  between  ma  and  p2  in  compressors.  To  find  B2,  p2 
was  plotted  versus  ma  at  constant  ng  values.  From  this  the  slope 
9ma/9p2  (B2)  was  obtained.  The  B1  expression  (9ma/9ng)  was  found 
by  regression  using  the  plotted  value  of  B2  in  the  Minitab  routine.  A2 
and  C2  were  then  found  using  plotted  values  to  find  the  slopes  9t2/9p2 
(A2)  and  9qc/9p2  (C2).  These  were  then  used  as  inputs  into  the 
relationships: 

B1  =  9ma/9ng  =  A1  =  9t2/9ng  =  Cl  =  9qc/9ng 

B2  9ma/9p2  A2  9t2/9p2  C2  9qc/9p2 

from  which  A1  and  Cl  were  obtained.  The  results  were  then  checked 
and  verified  through  the  use  of  the  data  regression  program.  At  steady 
state  qc  =  qhpt  is  required,  thus  the  relationship  for  Cl  (9qc/9ng)  =  14 
(9qhpt/9ng).  So.  with  14  determined,  the  value  of  II  (9qhpt/9ma)  can 
be  found  as  II  =  I4/B1  and  12  (9qhpt/9t2)  is  I4/A1.  With  these  three. 
II,  12,  and  14,  as  inputs,  13  (9qhpt/9p4)  and  15  (9qhpt/9mf)  were 
obtained  through  the  use  of  the  data  regression  routine.  The 
numerical  values  can  be  found  in  Appendix  C. 

The  next  step  involves  fine  tuning  the  computer  simulation  model 
and  results  in  the  setting  of  the  remaining  coefficients  of  the  A  and  B 


matrices. 


In  fine  tuning,  the  actual  dynamic  data  was  used  in  conjunction  with 
a  simulation  approach  to  select  final  values  for  the  model  constants. 

The  simulation  approach  used  the  simple  linear  state  space 
equations  for  dng,  dns,  and  de.  These  described  a  first  order 
differential  equation  set,  with  time  varying  A  and  B  matrices,  to  mimic 
the  actual  nonlinearity  of  the  system  response.  This  was  accomplished 
by  approximating  the  true  nonlinear  relationships  with  a  series  of 
small  linear  steps  as  depicted  in  figure  7. 

The  dynamic  simulation  computer  model  was  developed  from  a 
routine  written  by  Herda  and  later  modified  by  Stammetti,  it  is  shown 
in  figure  8.  The  major  change  to  the  previous  program  was  in  the  A 
and  B  matrix  equations.  In  the  past,  the  equations  were  complex 
and  were  a  combination  of  exponentials  and  second  order  equations. 
In  the  present  model  they  were  simple,  first  order  linear  functions  of 
the  states.  Another  change  which  occurred  in  the  fine  tuning  process 
was  the  adjusting  of  the  dynomometer  shaft  inertia.  This  change 
was  for  shaping  the  response  and  did  not  affect  the  steady  state 
response  of  the  model. 

A.  Computer  Simulation 

The  simulation  program  was  built  using  four  sections.  These  were 
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Time 

Figure  7.  Sequential  Linearization. 


Input: 


A(t),  B(t) 


Compute: 


x  =  A(t)  x  +  B(t)  u 


Figure  8.  Simulation  Process. 


the  dynamic  data  input,  the  initial  condition,  the  dynamic,  and  the 
derivative  sections.  The  actual  dynamic  data  was  taken  on  a  four 
channel  strip  chart  recorder  and  depicts  the  real  turbine  response  of 
Ng,  Ns,  Qd,  and  Mf  as  a  function  of  time.  These  dynamic  runs  were 
compiled  in  Appendix  E  and  provided  the  source  and  data  for  the 
simulation.  The  initial  condition  section  established  the  initial 
conditions  for  each  particular  run.  The  dynamic  section  formulated 
the  display  of  the  dynamic  data  and  was  followed  by  the  derivative 
section  which  formulated  the  simulated  response.  The  outputs  of 
these  final  two  sections  were  then  plotted  for  comparison.  These 
plotted  remits  were  used  to  assist  in  a  grid  search  on  the  still  doubtful 
component  coefficients  to  determine  the  final  values  of  the  A  and  B 
matrices.  The  procedure  was  to  guess  values  of  the  coefficients,  then 
simulate  the  model  response.  If  the  model  and  the  data  agreed,  then 
the  search  was  over. 

B.  The  Grid  Search 

A  grid  search  technique  was  applied  in  two  and  three  parameter 
searches  to  find  the  optimal  values  for  the  remaining  coefficients.  It 
was  formed  using  a  the  spread  sheet  formatted  program  called  Excell, 
by  Microsoft  software.  An  example  of  this  is  found  in  Appendix  C.  The 
top  portion  provided  the  inputs  for  the  coefficients  used  in  the 
calculation  of  the  A  and  B  matrices.  The  middle  portion  calculated  the 
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intermediate  Y  variables.  These  were  then  combined  with  the 
previously  mentioned  A  thru  K  coefficients  to  calculate  the  A  and  B 
matrix  values,  three  of  which  were  ng  dependent. 

We  also  knew  that  the  steady  state  required  x  =  0.  This  fact  was 
used  to  generate  additional  equations  which  must  be  satisfied  by  any 
coefficient  set  to  guide  the  grid  search.  The  last  section  of  the  spread 
sheet  provided  for  the  calculation  of  the  values  of  dng  and  dns  for  the 
simulation  runs  3  and  9  as  well  as  the  peak  dns  value  for  run  three 
(figure  9).  Run  three,  shown  in  figure  9,  developed  into  the  most 
difficult  run  for  determining  the  best  coefficients  for  Ns  simulation 
and  is  the  reason  for  the  addition  of  the  peak  dns  value.  The  values 
for  the  G,  F2,  and  F3  coefficients  were  the  focus  of  the  grid  search  as 
we  had  least  confidence  in  them.  These  are  the  important  coefficients 
in  determining  the  value  of  dns  and  in  turn,  the  ns  simulation  values. 
The  G1  -  G5  values  were  calculated  in  the  same  manner  as  described 
eariler  for  the  11-15  values.  They  represent  the  partial  derivatives  of 
the  equation  for  t4  (eqn.  10)  and  were  determined  to  have  little  effect 
on  the  A21.  A22  and  A33  results.  Their  elimination  as  a  search 
parameter  reduced  the  search  to  a  two  parameter  one.  Various  values 
for  F2  and  F3  were  tried,  and  the  resulting  A  and  B  matrix  coefficients 
calculated  were  inserted  in  the  computer  simulation.  The  results 
were  plotted  on  the  NPS  mainframe  computer  system,  using  the 
easyplot  routine,  to  check  for  accuracy  of  simulation.  This  process  was 
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Figure  9.  Dynamic  Data  Runs. 


repeated  until  the  optimal  F2  and  F3  values  were  obtained. 
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The  validation  of  the  model  was  a  two  step  process  which  carried 
over  from  the  fine  tuning.  After  the  fine  tuning  was  complete,  and 
the  model  resulted  in  the  proper  shape  for  runs  3  and  9,  and  with  the 
steady  state  response  within  the  required  tolerances,  the  first  step 
was  complete.  The  next  step  was  to  take  the  model  and  apply  it 
elsewhere  in  the  operating  range  of  the  turbine  to  check  for  accuracy. 
If  the  proper  response  in  shape  and  prescribed  tolerance  criteria  (+/ 
-  10  percent  in  the  dynamic  and  steady  response)  were  still  met, 
then  the  validation  was  completed.  This  process  was  only  performed 
around  the  borders  of  the  operating  range  (figure  9).  If  successful, 
then  the  model  should  be  applicable  anywhere  in  the  turbine 
operating  range. 

The  first  run  modeled  was  the  one  for  run  nine  which  models  the 
Ns  transient  response  from  500  to  2500  rpms  while  holding  the  value 
for  Ng,  as  load  varies,  constant  at  32000  rpms  (+/-  5%)  for  the 
conditions  shown  in  run  nine  in  Appendix  E..  Figure  10  is  the 
graphical  representation  of  the  modeled  run  and  the  dynamic  data  run 
for  this  transient.  The  criteria  has  been  met  for  run  9,  so  the  model 
was  applied  to  the  validation  run  one.  Run  one,  simulates  the 
transient  for  Ns  from  500  to  2215  rpms  (the  turbine  will  not  operate 
above  2215  rpm  at  low  Ng)  with  the  value  for  Ng  constant  at  20000 
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SIMULATION  RUN  9  NS 


Figure  10.  Simulation-Run  Nine. 
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rpms  (+/-  5%);  the  results  are  shown  in  figure  11.  Again,  the 
accuracy  criteria  were  met.  This  covers  the  upper  and  lower  bounds 
of  the  operating  range. 

Next,  run  three  was  used  in  conjunction  with  run  nine  in  the  fine 
tuning  portion.  This  run  is  a  vertical  run  on  the  left  side  of  the 
operating  range  and  is  an  Ng  transient  from  20000  to  32000  rpms 
with  Ns  being  held  constant  at  500  rpms  to  conform  with  the  actual 
dynamic  data  run.  This  run  is  shown  in  figure  12.  Note  that  the 
accuracy  criteria  are  easily  met  in  Ng  but  apparently  not  satisfied  in 
Ns.  However,  when  the  accuracy  criteria  is  applied  to  the  normalized 
Ns.  the  transient  is  acceptable.  To  validate  this  run  the  model  is  next 
applied  to  the  right  hand  side  of  the  operating  range.  This  is  run 
seven  and.  again,  is  a  Ng  transient  from  20000  to  32000  rpms  with  Ns 
slightly  increasing  from  the  lower  constraint  of  2215  rpm  to  2500 
rpms  which  is  attainable  at  the  upper  limits  of  the  operating  range. 
This  run  is  presented  in  figure  13.  A  validation  run  was  attempted 
vertically  in  the  center  of  the  operating  range.  This  is  run  five  and 
it's  response  was  of  the  same  result  shown  in  figure  1 1  for  run  seven 
and  therefore  was  not  shown  as  a  separate  figure. 
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SIMULATION  RUN  3  NG 
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Figure  12.  Simulation-Run  Three. 
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VALIDATION  RUN  7  NG 
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Figure  13.  Validation-Run  Seven 


VI. 


A  modeling  method  has  been  developed  for  real  time  simulation  of  a 
marine  gas  turbine.  The  results  of  this  modeling  method  should  be 
used  with  the  LQR  controller  design  theory  to  develop  a  new 
controller  for  the  Boeing  test  bed  at  NPS.  This  controller  should  then 
be  tested  to  ascertain  it's  ability  to  be  more  fuel  efficient  than  the 
present  system. 

This  model  was  simple  in  that  it  used  matrix  elements  for  A  and  B 
which  were  only  linear  state  functions.  The  dependency  was  only  on 
the  Ng  state.  This  simplicity  will  provide  fast  run  times  on  a  small 
computer  for  model  reference  control  and  diagnostics. 

Accurate  simulation  performance  has  been  demonstrated  in  the 
compressor/high  pressure  turbine  section  and  the  resulting  steady 
states  have  been  achieved  to  +/-  3%.  An  exception  to  this  accuracy  is 
the  constant  Ns  simulation  for  the  vertical  runs  .  This  is  only  an 
indication  that  the  final  values  for  F2,  F3,  and  possibly  G1  are  not  fine 
tuned  enough  to  provide  for  more  accurate  constant  Ns  response. 
Investigation  was  conducted  into  this  irregularity  to  determine  the 
cause  of  the  steady  state  results  being  to  high.  The  fact  that  the  steady 
state  response  for  a  constant  Ns  is  higher  in  runs  five  and  seven  was 
probed  and  a  number  of  solutions  were  looked  at.  One  such  solution 
was  to  check  for  the  possibility  that  the  oscillating  torque  for  runs  five 
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and  seven  shown  in  Appendix  E  was  causing  the  Ns  value  to  continue 
to  rise  and  settle  out  at  a  higher  steady  state.  These  oscillations  do 
not  appear  to  be  the  cause  nor  have  any  affect  on  the  steady  state 
response  of  Ns.  This  was  verified  by  modifying  the  computer 
simulation  routine  so  as  to  remove  the  oscillation.  Another  solution 
was  to  adjust  the  A21  values  until  a  proper  steady  state  value  for  ns  was 
reached.  At  the  time  of  publication  this  has  not  provided  any 
satisfactory  results.  Another  possibility  is  that  the  A21  coefficients 
may  have  to  change  the  sign  to  produce  the  proper  response.  The  last 
area  to  examine  as  a  possible  solution  to  this  problem  is  to  input  a 
varying  inertia  for  the  water  brake  dynamometer.  In  the  past  this 
inertia  has  been  assumed  to  be  a  constant  value  when,  in  fact,  it 
actually  changes  with  the  addition  and  subtraction  of  water. 

The  last  recommendations  for  this  project  encompass  the  entire 
project.  First,  a  digital  data  acquisition  system  should  be  used  to  get 
the  steady  state  and  dynamic  data.  This  will  make  the  fine  tuning 
process  with  the  data  regression  program  more  accurate.  Finally,  one 
computer  system  should  be  used  for  the  whole  process,  from  the 
regression  and  grid  search  to  the  computer  simulation  (vice  the  two 
system  used  in  this  project).  In  this  way,  results  could  be  achieved 
and  correlated  in  a  more  efficient  manner. 

The  model  was  robust  in  that,  with  the  exception  of  the  constant 
Ns  in  the  vertical  runs,  the  model  was  applicable  to  all  areas  of  the 
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operating  range.  This  was  supported  by  the  excellent  response  in  Ng 
throughout  the  whole  operating  range.  This  shows  that  the 

coefficients  for  the  All,  A12.  and  Al3  were  accurate  and  the 
procedure  to  find  them  was  sound. 

Once  the  proper  values  for  A21,  A22,  and  A23  are  found  and  this 
modeling  technique  has  been  successfully  demonstrated  on  the  Boeing 
engine  this  method  will  have  real  time  application  to  the  future 
controller  design  for  the  General  Electric  LM-  2500  and  the  follow  on 
marine  gas  turbine  controllers  designed  for  the  future. 
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APPENDIX  A:  COMPONENT  MODELING  EQUATIONS 

APPENDIX  B:  EQUATION  SET  REDUCTION  METHODOLOGY 

APPENDIX  C:  PARAMETER  SELECTION  METHODOLOGY 

APPENDIX  D:  COMPUTER  SIMULATION  CODES 

APPENDIX  E:  SOURCE  DATA  RUNS 


AEEBHPIS  A 


COMPONENT  MODELING  EQUATIONS  : 


Compressor 


dt2  =  Al  dng  +  A2  dp2 
dma  =  Bl  dng  +  B2  dp2 
dqc  =  Cl  dng  +  C2  dp2 


1 

2 

3 


Gas  Generator  Shaft 


JGB  dng  =  dqhpt-dqc 


Dynamometer  Shaft 


JDB  dns  =  dqfpt-dqd 
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dp4  =  D1  dng  +  D2  dns  .  6  * 


dqfpt  =  Fl  dns+F2  dt4+F3  dmaf  ...  7 


High  Pressure  Turbine 


10 
11  * 


dt4  =  Gl  dma  +  G2  dt2  +  G3  dp4  +  G4  dng  +  G5  de 

dp2  =  Hi  dng  +  H2  dns  . 

dqhpt  =  U  dma  +  12  dt2  +  13  dp4  +  14  dng  +  15  de 


12 


Note:  1.  Lower  case  variables  represent  normalized  variables. 
2.  The  d  _  represents  the  perturbation  of  the  variable. 

*  Decoupling  equations. 
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EQUATION  SET  REDUCTION  METHODOLOGY: 


The  equations  in  appendix  A  are  reduced  into  the  state  space 
components  of  the  A  and  B  matrices. 


Substitute  eqn  1 1 : 

dp2  =  Hi  dng  +  H2  dns  .  1 1 

into  eqns  1,  2,  and  3  for  dp2. 

dt2  =  (Al  +  A2Hl)  dng  +  (A2H2)  dns . 1 

dma  =  (Bl  +  B2H1)  dng  +  (B2H2)  dns . 2 

dqc  =  (Cl  +  C2H1)  dng  +  (C2H2)  dns . 3 

In  eqn  12: 

dqhpt  =  II  dma  +  12  dt2  +  13  dp4  +  14  dng  +  15  de  .  12 

substitute  eqns  1,  2,  and  6  for  dt2,  dma,  and  dp4. 
dqhpt  =  11  (Bl  +  B2H1)  dng  +  II  (B2H2)  dns  + 

12  (Al  +  A2H1)  dng  +  12  (A2H2)  dns  + 

13  (D  l )  dng  +  13  (D2)  dns  + 

14  dng  + 15  de  .  12 
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Collect  terms: 


dqhpt  =  (II  (Bl  +  B2H1)  +  12  (Al  +  A2H1)  +  13  (Dl)  +  I4  )  dng  + 
(U  (B2H2)  +  12  (A2H2)  +  13  (D2)  )  dns  +  15  de 

Let  Y1  =  II  (Bl  +  B2H1)  +  12  (Al  +  A2H1)  +  13  (Dl)  +  14 
Y2  =  U  (B2H2)  +  12  (A2H2)  +  13  (D2) 

Then:  dqhpt  =  Y1  dng  +  Y2  dns  +  15  de  .  12 

Take  eqn  4  and  substitute  in  new  eqn  12  for  dqhpt  and  eqn  3  for 
dqc. 

JGB  dng  =  Y1  dng  +  Y2  dns  +  15  de  - 

(Cl  +  C2H1)  dng  -  (C2H2)  dns 
Collect  terms: 

dng  =  (  (Y1  -  (Cl  +  C2H1))  dng  +  (Y2  -  (C2H2))  dns  +  15  de  )/JGB 

Let:  Z1  =  (Y1  -  (Cl  +  C2H1) 

Z2  =  Y2  -  (C2H2) 

Z3  =  15 

Then: 

dng  =  (Z1  dng  +  Z2  dns  +  Z3  de)/JGB  . R1 
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In  eqn  10: 

dt4  =  Gl  dma  +  G2  dt2  +  G3  dp4  +  G4  dng  +  G5  de  . . .  10 

substitute  eqns  1,  2,  and  6  for  dt2,  dma,  and  dp4. 
dt4  =  Gl  (B1  +  B2H1)  dng  +  Gl  (B2H2)  dns  + 

G2  (Al  +  A2Hl)  dng  +  G2  (A2H2)  dns  + 

G3  (Dl)  dng  +  G3  (D2)  dns  + 

G4  dng  +  G5  de  .  10 

Collect  terms: 

dt4  =  (Gl  (Bl  +  B2H1)  +  G2  (Al  +  A2Hl)  +  G3  (Dl)  +  G4  )  dng  + 
(Gl  (B2H2)  +  G2  (A2H2)  +  G3  (D2)  )  dns  +  G5  de 

Let  Y3  =  Gl  (Bl  +  B2H1)  +  G2  (Al  +  A2Hl)  +  G3  (Dl)  +  G4 
Y4  =  Gl  (B2H2)  +  G2  (A2H2)  +  G3  (D2) 

Then:  dt4  =  Y3  dng  +  Y4  dns  +  G5  de  .  10 

In  the  continuity  eqn  (9)  substitute  eqn  2  for  dma. 

dmaf  =  Kl  dmf  +  K2(Bl  +  B2H1)  dng  +  K2(B2H2)  dns  . 9 

Substitute  new  eqn  10  for  dt4  and  new  eqn  9  for  dmaf  into: 
eqn  7:  dqfpt  =  Fl  dns  +  F2  dt4  +  F3  dmaf 
This  results  in: 
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dqfpt  =  Fl  dns  +  F2V3  dng  +  F2Y4  dns  +  F2G5  de  + 

F3K1  dmf  +  F3K2(B1  +  B2H1)  dng  +  F3K2(B2H2)  dns 
Collect  terms: 

dqfpt  =  (F2Y3  +  F3K2(Bl  +  B2Hl))  dng  + 

(Fl  +  F3K2(B2H2)  +  F2Y4)  dns  +  F2G5  de  +  F3K1  dmf 

Let:  Z4  =  F2Y3  +  F3K2(Bl  +  B2H1) 

Z5  =  Fl  +  F3K2(B2H2)  +  F2Y4 
Z6  =  F2G5 
Z7  =  F3K1 
Then: 

dqfpt  =  Z4  dng  +  Z5  dns  +  Z6  de  +  Z7  dmf  . 7 

Substitute  new  eqn  7  for  qfpt  in  eqn  5. 


JDB  dns  =  dqfpt  -  dqd  . 5 

dns  =  (Z4  dng  +  Z5  dns  +  Z6  de  +  Z 7  dmf  -  dqd) /JDB  . R2 


Finally:  Eqn  8  is  R3.  Rl,  R2,  and  R3  make  up  the  state  space 
equation  set. 
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The  state  space  equation  set  is  the  following: 


dng  =  21  dng  +  22  dns  +  25  de  . R1 

JGB  JGB  JGB 

dns  =  2A  dng  +  25  dns  +  Z6  de  +  Z7  dmf  -  X  dqd  . R2 

JDB  JDB  JDB  JDB  JDB 

de  =  -  1  de  +  1  dmf  . R3 

T  T 


The  resulting  A  and  B  matrices  are: 


All 

=  Zi 

A12  =  Z2 

A13 

=  25 

JGB 

JGB 

JGB 

A21 

=  Z4 

A22  =  Z5 

A23 

=  Z6 

JDB 

JDB 

JDB 

A31 

=  O 

A32  =  O 

A3  3 

= 

T 


and 

Bll  =  O  B12  =  O 


B2 1  =  Z7  B22  =  -1 
JDB  JDB 

B31  =  _L  B32  =  O 
T 


Afpgypg  £ 

PARAMETER  SELECTION  METHODOLOGY: 


The  parameter  selections  were  accomplished  through  three  means: 

1 .  Plotting  actual  slopes  from  data.  This  was  shown 
in  chapter  III  section  F. 

2.  Data  reduction  using  mainframe  program 
Mini-Tab. 

3.  Grid  searching  using  coefficients  from  1.  and  2. 
and  known  relationships  which  must  be  satisfied. 
The  spread  sheet  is  from  the  Micro-soft  program 

Excell  and  was  used  on  an  Apple  Macintosh 
computer. 

The  data  reduction  program  Minitab  was  used  for  all  of  the 
regressions  performed.  The  input  for  the  regression  program  was  the 
steady  state  turbine  data  shown  in  table  1  A.  This  routine  was  used  to 
determine  the  order  of  the  independent  variables  in  the  equations  for 
Ma,  T2,  and  Qc.  as  well  as  to  regress  the  expressions  for  P2,  P4,  T4, 
Qhpt  and  Maf.  To  determine  the  order  of  a  independent  variable  a 
three  step  process  was  followed.  First,  the  dependent  variable  was 
regressed  using  the  data  for  the  independent  variables.  In  doing  this  a 
residual  was  calculated  for  each  set  of  the  independent  variables.  This 
was  then  subtracted  from  the  regressed  value  of  the  dependent 
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variable.  The  resulting  difference  was  plotted  versus  each  of  the 
independent  variables.  In  the  plot  for  the  regression  of  T2,  shown  in 
figure  C-la,  there  is  a  definite  parabolic  shape  which  is  an  indication 
that  the  independent  variable  in  question  (Ng)  is  of  a  higher  order.  To 

verify  this,  T2  was  regressed  again,  this  time  using  NgA2,  Ng  and  P2. 
The  plotting  was  repeated  with  the  resultant  plot  shown  in  figure  C- 
lb.  This  time  the  parabolic  shape  was  not  in  evidence  and  the 
predictor  percentage,  which  indicates  the  accuracy  of  the  regression 
based  on  the  data,  increased  from  99.6  to  99.9  percent  indicating  an 
excellent  regression  for  the  t2  equation.  The  resulting  partials  for  all 
of  the  component  expressions  are  shown  in  figure  C-2  lines  3  through 
21. 

The  spreadsheet  grid,  figure  C-2,  was  used  in  conjunction  with  the 
simulation  program  to  search  out  the  best  values  for  the  rest  of  the 
coefficients  in  the  equations  in  appendix  A.  These  were  input  for  Gl, 
F2,  and  F3  in  lines  3-21  and  then  used  in  lines  25-32  to  calculate  the 
values  of  the  coefficients  of  the  A  and  B  matrices.  Lines  34-38  show 
the  proper  signs  for  these  matrix  coefficients  as  well  as  the  dns  and 
dng  relationships  which  must  be  satisfied  during  the  search  process. 
Once  these  were  satisfied,  plots  were  generated  and  checked  for  the 
accuracy  of  the  simulation. 
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Figure  C-l.  Regression  Plots 
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*  * 


Figure  C-2.  Sample  Grid  Search  Spread  Sheet. 
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COMPUTER  SIMULATION  CODES: 


The  computer  codes  were  written  using  the  Dynamic  Simulation 
Language  (DSL). 

Simulation  Run  1 :  DSL  code  for  dynamic  Ns  run  at  constant  Ng 

(20000  rpm). 

Simulation  Run  3:  DSL  code  for  dynamic  Ng  run  at  constant  Ns 

(500  rpm). 

Simulation  Run  5:  DSL  code  for  dynamic  Ng  run  at  constant  Ns 

(1500  rpm). 

Simulation  Run  7  :  DSL  code  for  dynamic  Ng  run  at  Ns 

(  2252-2500  rpm). 

Simulation  Run  9:  DSL  code  for  dynamic  Ns  run  at  constant  Ng 

(32000  rpm). 
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*  * 

*  * 


*  BOEING  MODEL  502 -6A  GAS  TURBINE 

* 


* 

* 


*  DYNAMIC  COMPUTER  SIMULATION 

* 


* 

* 


*  MODIFIED  BY 

*  S.D.  METZ 

*  FROM  ROUTINES  BY 

*  V. J. HERDA  AND  V.  A.  STAMMETTI 


* 

* 

* 

* 

* 


*  THIS  PROCRAM  SIMULATES  THE  DYNAMIC  RESPONSE  OF  THE  NPS  * 

*  BOEINC  GAS  TURBINE  TEST  FACILITY  USING  A  MULTILPLE  * 

*  LINEARIZATION  TECHNIQUE.  * 

*  * 


*********************************************************************** 


C 

C 

P ARAM  JG=0. 009525,  JD=0.  3000,  FI=3. 14159,  TC=1. 0 


THE  FOLLOWING  VALUES  LISTED  ON  THE  FUNCTION  CARD  ARE  FOR  FUEL  FLOW, 
GAS  GENERATOR  SPEED,  TORQUE  AND  DYNO  SPEED  AS  A  FUNCTION  OF  TIME. 
THESE  VALUES  WERE  OBTAINED  FROM  STRIP  CHART  RECORDS  AND  ARE  ENTERED 
IN  THE  FORM  ( E. G.  FUEL  FLOW)  . . . TIME( SEC) ,  FUEL  FLOW . 


C 

C  THIS  SET  IS  FOR  EXPERIMENTAL  RUN  #  1. 

C 

C 

AFGEN  NGDATA=  0. 0, 20230. 0,  5.0,20230.0,  10.0,20230.0,  15.0,20230.0,  ... 

20.0,20230.0,  25.0,20230.0,  50.0,20230.0 
AFGEN  NSDATA=  0. 0, 505. 0,  1.0,505.0,  2.0,536.7,  3.0,568.33,  ... 

4.0,600.0,  5.0,631.7,  6.0,695.0,  7.0,790.0,  8.0,885.0,  ... 
9.0,1043.3,  10.0,1170.0,  11.0,1328.3,  12.0,1486.7,  ... 
13.0,1708.3,  14.0,1835.0,  15-0,1930.0,  16.0,1993.3,  ... 

17.0,2056.7,  18.0,2088.3,  19.0,2120.0,  20.0,21 51.7,  ... 

21.0,2167.5,  22.0,2177.0,  23.0,2183.3,  24.0,2215.0,  ... 

25.0,2215.0,  50.0,2215.0 

AFGEN  MFDATA=  0. 0,78. 39,  5.0,78.39,  10.0,78.39,  12.0,79.89,  ... 
15.0,80.26,  17.0,80.26,  19.0,81.38,  20.0,82.50,  ... 

21.0,82.13,  25.0,82.13,  50.0,82.13 


AFGEN  QDDATA=  0. 0, 125. 0,  1.0,125.0,  2.0,118.3,  3.0,115.0, 
4.0,111.7,  5.0,105.0,  6.0,98.33,  7.0,85.0,  ... 
8.0,71.67,  9.0,58.33,  10.0,45,0,  11.0,31.67,  ... 
12.0,11.67,  12.5,1.670,  13.0,5.0,  15.0,5.0,  ... 
25.0,5.0,  50.0,5.0 


INITIAL 

*  ESTABLISH  INITIAL  CONDITIONS. 

* 

NGI  =  20230.  0 
NSI=505  0 
MFI  =  78.  39 
QDI  =  125. 0 


*  SET  INITIAL  STATE  PERTURBATION  TO  ZERO 


DNG  =  0.0 
DNS  =  0.  0 
DE  =  0.  0 
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NGB  =  26000. 0 

NSB  =  1500. 0 

MFB  =  122. 5 

QDB  =  220. 0 

QCB  =  80. 0 

JGB  =  (NGB/QCB)*JG 

JOB  =  ( NSB/QDB ) *  JD 

DYNAMIC 

* 

*  DATA  CURVE  FORMULATION 

* 

NGD  =  AFGEN( NGDATA , T I ME ) 

NSD  =  AFGEN(  NSDATA , T I ME ) 

MFD  =  AFGEN( MFDATA , T I ME ) 

QDD  =  AFG£N( QDDATA , T I ME ) 

MFDLY  =  TRANSP(  100,MFI,0.  70, MFD) 

* 

*  STATE  SPACE  LINEAR  MODEL  FORMULATION 

* 

All  =  (20. 282-32. 560*( NGL/NGB ) )/(  JGB) 

* 

A12  *  ( -0. 4480 )/( JGB) 

* 

A13  =  ( 13. 624)/( JGB) 

* 

A21  =  (2. 1131+.91786*{ NGL/NGB) )/( JDB) 

* 

A22  =  (-. 0998-0. 445*( NGL/NGB) )/( JDB) 

* 

A23  =  (-0. 1885)/( JDB) 

* 

A3 1  =  0. 0 

* 

A32  =  0. 0 

* 

A33  =  -1.0/TC 

* 

Bll  =  0.0 

* 

B12  =  0. 0 

* 

B21  =  .  035250/( JDB) 

* 

B22  =  -1.  0/( JDB) 

* 

B31  =  1. 0/TC 

* 

B32  =  0. 0 


DERIVATIVE 

NOSORT 

* 

*  COMPUTE  INPUT  TO  THE  NONLINEAR  MODEL,  MF( T ) ,  WW(T). 

*  RUN  #9 

DMF  =( MFDLY  -  MF I ) /MFB 
DQD  =( QDD  -  QDI)/QDB 

DNGDOT  =  All *DNG  +  A12*DNS  +  A13*DE  +  B11*DMF 

DNSDOT  =  A2 1 *DNG  +  A22*DNS  +  A23*DE  +  B2 1 *DMF  ♦  B22*DQD 

DEDOT  =  A33*DE  *  B31*DMF 


*  DYNAMIC  EQUATIONS  FOB  LINEAR  MODEL. 
DNG=INTGRL( 0. 0, DNGDOT) 

DNS= INTGPL( 0.  0, DNSDOT) 

DE  =INTGRL( 0. 0, DEDOT) 

NCL  =  NGI  +  DNG*NGB 
NSL  =  NS  I  +  DNS  *NSB 

SORT 
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* 

*  THE  STATEMENTS  IN  THE  PREVIOUS  (DERIVATIVE)  SECTION  YIELD  VALUES 

*  OF  ’NO’,  AND  'NS'  AS  CALCULATED  BY  THE  NONLINEAR  AND  STATE  SPACE 

*  MODELS.  THE  STATEMENTS  BELOW  COMPUTE  THE  VALUES  OF  'NO',  AND  'NS' 

*  AS  RECORDED  FROM  CAS  TURBINE  TEST  DATA. 

* 

TERMINAL 
METHOD  RK5 

CONTROL  FINTIM=50. 0,DELT=0. 001 
PRINT  0.  S^GD.NGL.NSD^SL/MFD/QDD 
END 
STOP 
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* 


* 


* 


* 


*  BOEING  MODEL  502-6A  GAS  TURBINE 


* 

* 


* 

* 


DYNAMIC  COMPUTER  SIMULATION  * 

* 


*  MODIFIED  BY 

*  S.  D.  METZ 

*  FROM  ROUTINES  BY 

*  V.  J. HERDA  AND  V. A.  STAMMETTI 


* 


* 


*  THIS  PROGRAM  SIMULATES  THE  DYNAMIC  RESPONSE  OF  THE  NPS  * 

*  BOEING  GAS  TURBINE  TEST  FACILITY  USING  A  MULTILPLE  * 

*  LINEARIZATION  TECHNIQUE.  * 

*  * 


*********************************************************************** 


C 

C 

PARAM  JC=0. 009525,  JD=0. 3000,  PI=3. 14159,  TC=1.  00 

* 

*  THE  FOLLOWING  VALUES  LISTED  ON  THE  FUNCTION  CARD  ARE  FOR  FUEL  FLOW, 

*  GAS  GENERATOR  SPEED,  TORQUE  AND  DYNO  SPEED  AS  A  FUNCTION  OF  TIME. 

*  THESE  VALUES  WERE  OBTAINED  FROM  STRIP  CHART  RECORDS  AND  ARE  ENTERED 

*  IN  THE  FORM  (E.G.  FUEL  FLOW)  .  . . TIME(  SEC) ,  FUEL  FLOW . 

* 


C 

C  THIS  SET  IS  FOR  EXPERIMENTAL  RUN  #  3. 

C 

c 

AFC-EN  NGDATA=  0.  0,20220.  0,  3.55,20220.0,  3.7,20315.7,  3.8,20411.4,  ... 

4.0,20698.4,  5.0,23186.1,  6.0,25482.4,  7.0,27204.6,  ... 
8.0,28926.9,  9.0,30649.1,  10.0,32180.0,  10.5,32371.4,  ... 
11.0,32180.0,  11.5,32084.3,  12.0,32180.0,  50.0,32180.0 
AFGEN  NSDATA=  0. 0, 506. 0,  2.0,505.7,  3.6,505.5,  4.0,506.5,  ... 

5.0,508.3,  6.0,509.7,  7.0,511.1,  8.0,512.5,  9.0,513.5,  ... 
10.0,514.4,  10.4,514.7,  11.0,514.4,  12.0,515.1,  ... 

13.0,515.8,  14.0,516.7,  16.3,527.2,  27.0,526.7,  ... 

18.0,515.8,  19.0,514.9,  20.0,513.0,  50.0,513.0 

AFGEN  MFDATA=  0. 0, 77. 5,  3.0,77.5,  3.3,79.6,  3.6,85.8,  ... 

4.0,91.99,  5.0,106.5,  6.0,118.9,  7.0,130.5,  8.0,141.3,  ... 
9.0,156.2,  10.0,172.7,  11.0,156.2,  12.0,160.3,  50.0,160.3 


0,  3.8,  128.  4,  4.0,131.  8,  ... 
0,274.1,  8.0,321.6,  9.0,375.8, 
11.0,430.0,  1 2.0,416.4,  ... 
15.0,416.4,  18.0,423.0,  ... 


INITIAL 

1  ESTABLISH  INITIAL  CONDITIONS. 

t 

NGI=20220. 0 
NSI=506. 0 
MFI  =  77. 5 
QDI  =  125. 0 


AFGEN  QDDATA=  0.0,125.0,  3.6,125. 

5.0,172.4,  6.0,226.7,  7. 
10.0,430.0,  10.5,443.6, 
13.0,416.4,  14.0,419.8, 
20.0,430.0,  50.0,430 


*  SET  INITIAL  STATE  PERTURBATION  TO  ZERO 

* 

DNG  =  0.0 
DNS  =  0.0 
DE  =0.0 
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M 


NGB  =  26000.0 

NSB  =  1500.0 

MFB  =  122. 5 

QDB  =  220. 0 

QCB  =  80. 0 

JGB  =  (  NGB/QCB) * JG 

JOB  =  ( NSB/QDB ) *  JD 

DYNAMIC 

* 

*  DATA  CURVE  FORMULATION 

• 

NGD  =  AFGEN( NGD ATA , T I ME ) 

NSD  =  AFGEN( NSDATA, TIME ) 

MFD  =  AFGEN ( MFDATA , T I ME ) 

QDD  =  AFGEN( QDDATA, TIME ) 

MFDLY  =  TRANSP(  100, MFI ,0. 70, MFD) 

* 

*  STATE  SPACE  LINEAR  MODEL  FORMULATION 


* 


* 


* 


* 

* 

* 

* 

* 


* 


* 


All  =  (20. 282-32. 560*( NGL/NGB) )/(  JGB) 
A12  =  ( -0. 448 )/( JGB) 

A13  =  ( 13. 624 ) /( JGB) 

A21  =  ( 2. 1131+. 91786*( NGL/NGB) )/< JDB) 
A22  =  (-. 0998-0. 445*( NGL/NGB ))/(  JDB) 
A23  =( -. 18851 )/( JDB) 

A3 1  =  0. 0 

A32  =  O. 0 

A3 3  =  -1. 0/TC 

Bll  =  0. 0 

B12  =  0. 0 

B21  =. 035250/( JDB) 

B 22  =  -1. 0/< JDB) 

B31  =  1. 0/TC 
B32  =  0. 0 


DERIVATIVE 

NOSORT 

* 

*  COMPUTE  INPUT  TO  THE  NONLINEAR  MODEL,  MF( T) ,  WW(T). 

*  RUN  #3 

DMF  =( MFDLY  -  MFI )/MFB 
DQD  =( QDD  -  QDI )/QDB 

DNGDOT  =  All *DHG  +  A12*DNS  +  A13*DE  +  Bll *DMF 
DNSDOT  =  A21*DNG  +  A22*DNS  +  A23*DE  +  B21*DMF 
DEDOT  =  A33  *DE  +  B3 1 *DMF 


*  DYNAMIC  EQUATIONS  FOR  LINEAR  MODEL. 
DNG= I NTGP L( 0. 0 , DNGDOT ) 

DNS= I HTGRL( 0. O.DNSDCT) 

CE  =INTGRL( 0. 0, DEDOT) 

NCL  =  NGI  ♦  DHG*N3B 
NSL  =  NS  I  *  DNS*NSB 

SORT 


B22*DQD 
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*  THE  STATEMENTS  IN  THE  PREVIOUS  (DERIVATIVE)  SECTION  YIELD  VALUES 

*  OF  'NG',  AND  'NS'  AS  CALCULATED  BY  THE  NONLINEAR  AND  STATE  SPACE 

*  MODELS.  THE  STATEMENTS  BELOW  COMPUTE  THE  VALUES  OF  'NG',  AND  'NS* 

*  AS  RECORDED  FROM  GAS  TURBINE  TEST  DATA. 

* 

TERMINAL 
METHOD  RK5 

CONTROL  FINTIM=50. 0,DELT=0.  001 

PRINT  0.50, NGD , NGL , NSD ( NSL , MFD ( QDD , DNSDOT , A2 1 , DNG ,  ... 

A22 , DNS, A23 , DE, B21 , DMF, B22 , DQD 

END 

STOP 
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*  * 


*  * 

*  BOEING  MODEL  502-6A  GAS  TURBINE  * 

*  * 

*  DYNAMIC  COMPUTER  SIMULATION  * 

*  * 


*  MODIFIED  BY  * 

*  S. D.  METZ  * 

*  FROM  ROUTINES  BY  * 

*  V.  J. HERDA  AND  V. A.  STAMMETTI  * 

*  * 


*  THIS  PROGRAM  SIMULATES  THE  DYNAMIC  RESPONSE  OF  THE  NPS  * 

*  BOEING  GAS  TURBINE  TEST  FACILITY  USING  A  MULTILPLE  * 

*  LINEARIZATION  TECHNIQUE.  * 

*  * 


C 

C 

PARAM  JG=0. 009525,  JD=0. 3000,  PI=3. 14159,  TC=1. 00 

* 

*  THE  FOLLOWING  VALUES  LISTED  ON  THE  FUNCTION  CARD  ARE  FOR  FUEL  FLOW, 

*  GAS  GENERATOR  SPEED,  TORQUE  AND  DYNO  SPEED  AS  A  FUNCTION  OF  TIME. 

*  THESE  VALUES  WERE  OBTAINED  FROM  STRIP  CHART  RECORDS  AND  ARE  ENTERED 

*  IN  THE  FORM  ( E. G.  FUEL  FLOW)  . . . TIME( SEC ) ,  FUEL  FLOW . 


C 

C  THIS  SET  IS  FOR  EXPERIMENTAL  RUN  #  5. 

C 

C 

AFGEN  NGDATA=  0. 0,20220. 0,  3.1,20220.0,  3.6,20412.1,  4.0,21372.6,  . 
5.0,23677.7,  6.0,25886.9,  7.0,27711.8,  8.0,29440.6,  ... 
9.0,30977.4,  10.0,32206.8,  10.2,32322.1,  10.8,32130.0,  ... 

11.5,32130.0,  12.1,32014.7,  13.4,32206.8,  14.5,32014.7,  ... 

15.6,32206.8,  16.8,32014.7,  17.0,32206.8,  50.0,32130.0 

* 

AFGEN  NSDATA=  0. 0, 1491. 6,  1.0,1490.4,  3.0,1469.6,  4.0,1492.5,  ... 

5.0,1501.0,  6.0,1508.3,  7.0,  1511.7,  7.8,1513.3,  ... 
9.6,1513.3,  10.1,1513.8,  12.0,1508.3,  13.0,1510.0,  ... 

50.  0,1510. 0 

* 

AFGEN  MFDATA=  0. 0, 83. 75,  3.2,83.75,  5.0,109.6,  6.0,122.2,  ... 
7.0,135.5,  8.0,150.2,  9.0,162.8,  9.8,172.4,  ... 

10.0,168.7,  11.0,165.0,  50.0,165.0 

* 

AFGEN  QDDATA=  0.  0,71.  0,  4.4,71.0,  5.0,78.7,  6.0,94.1,  ... 

7.0,107.5,  8.0,120.9,  9.0,134.4,  10.2,149.8,  ... 
10.9,147.9,  11.6,147.9,  12.4,143.2,  13.0,145.9,  ... 
13.5,144.0,  14.0,145.9,  15.0,144.0,  50.0,144.0 


INITIAL 

ESTABLISH  INITIAL  CONDITIONS. 

NGI=20220. 0 
NSI =1491 . 6 
MFI  =  83. 75 
QDI  =71.0 


SET  INITIAL  STATE  FERTURBATION  TO  ZERO 


DNG  =  0.0 
DNS  =  0. 0 
DE  =  0.  0 
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NGB  =  26000. 0 

NSB  =  1500.0 

MFB  =  122. 5 

QDB  =  220. 0 

QCB  =  80. 0 

JGB  =  ( NGB/QCB ) * JG 

JDB  =  <NSB/QDB)*JD 

DYNAMIC 

* 

* 

* 

DATA  CURVE  FORMULATION 

ft 

NGD  =  AFGEN( NGDATA , T I ME ) 

NSD  =  AFGEN( NSDATA, TIME ) 

MFD  =  AFGEN( MFDATA, TIME ) 

QDD  =  AFGEN( QDDATA , T I ME ) 

• 

* 

MFDLY  =  TRANSPf 100, MF I „  0.  70, MFD) 

* 

* 

STATE  SPACE  LINEAR  MODEL  FORMULATION 

* 

All  =  ( 20. 282-32. 560*(  NGL/NGB) )/< JGB) 

* 

A12  =  ( -0. 4480 ) /(  JGB) 

* 

A13  =  ( 13. 624)/( JGB) 

* 

A21  =  ( 2. 1 13 1+ . 91786*( NGL/NGB ) ) /( JDB ) 

* 

A22  =  (-. 0998-0. 445*( NGL/NGB) )/(  JDB) 

* 

A23  =( -0. 18851 )/( JDB) 

* 

A31  =  0. 0 

* 

A32  =  0. 0 

* 

A33  =  -1.  0/TC 

* 

Bll  =  0. 0 

* 

B12  =  0. 0 

ft 

* 

B21  =. 035250/( JDB) 

* 

B22  =  -1. 0/( JDB) 

• 

* 

B3 1  =  1 .  0/TC 

★ 

B32  =  0. 0 

* 

DERIVATIVE 

NOSORT 

* 

* 

COMPUTE  INFUT  TO  THE  NONLINEAR  MODEL,  MF( T ) 

,  WW( T ) . 

* 

RUN  #9 

DMF  =( MFDLY  -  MF I ) /MFB 

DQD  =( QDD  -  QDI)/QDB 

DNGDOT  =  All *DNG  +  A12*DNS  +  A13*DE  + 
DNS DOT  =  A2 1 *DNG  +  A22*DNS  +  A23*DE  + 

Bll *DMF 

B21*DMF  *  E22*DQD 

* 

DEDOT  =  A33*DE  +  B31*DMF 

* 

* 

DYNAMIC  EQUATIONS  FOR  LINEAR  MODEL. 

DNG=INTC.RL(  0  0,  DNGDOT) 

DNS=INTGRL( 0. O.DNSDOT) 

DE  =INTGR'( 0. 0, DEDOT) 

NGL  =  NG I  +  DUG* NGB 

NSL  =  NS I  *  DHS*NSB 

CAD 

*T’ 

f 
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*  THE  STATEMENTS  IN  THE  PREVIOUS  (DERIVATIVE)  SECTION  YIELD  VALUES 

*  OF  'NG',  AND  ’NS'  AS  CALCULATED  BY  THE  NONLINEAR  AND  STATE  SPACE 

*  MODELS.  THE  STATEMENTS  BELOW  COMPUTE  THE  VALUES  OF  'NG' ,  AND  'NS' 

*  AS  RECORDED  FROM  GAS  TURBINE  TEST  DATA. 

* 

TERMINAL 
METHOD  RK5 

CONTROL  FINTIM=50. 0,DELT=0.  001 

PRINT  0.5, NGD , NGL , NSD , NSL , MFD , QDD , DNSDOT , A2 1 , DNG , A2 2 , DNS ,  ... 

A23 , DE, B2 1 , DMF, B22 , DQD 

END 

STOP 
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* 


* 


* 


* 


*  BOEING  MODEL  502-6A  GAS  TURBINE 

* 


* 

* 


*  DYNAMIC  COMPUTER  SIMULATION 

* 


* 

* 


*  MODIFIED  BY 

*  S.D.  METZ 

*  FROM  ROUTINES  BY 

*  V. J. HERDA  AND  V. A.  STAMMETTI 


* 


* 

* 

* 


*  THIS  PROGRAM  SIMULATES  THE  DYNAMIC  RESPONSE  OF  THE  NPS  * 

*  BOEING  GAS  TURBINE  TEST  FACILITY  USING  A  MULTILPLE  * 

*  LINEARIZATION  TECHNIQUE.  * 

*  * 


C 

C 

PARAM  JG=0. 009525,  JD=0.3000,  PI=3. 14159,  TC=1. 00 

« 

*  THE  FOLLOWING  VALUES  LISTED  ON  THE  FUNCTION  CARD  ARE  FOR  FUEL  FLOW, 

*  GAS  GENERATOR  SPEED,  TORQUE  AND  DYNO  SPEED  AS  A  FUNCTION  OF  TIME. 

*  THESE  VALUES  WERE  OBTAINED  FROM  STRIP  CHART  RECORDS  AND  ARE  ENTERED 

*  IN  THE  FORM  ( E.  G.  FUEL  FLOW)  . . , TIME( SEC ) ,  FUEL  FLOW . 


C 

C  THIS  SET  IS  FOR  EXPERIMENTAL  RUN  #  7. 

C 

C 

AFGEN  NGDATA=  0. 0, 20220. 0,  4.0,20220.0,  4.5,20759.7,  5.0,21954.7,  ... 

6.0,24460.3,  7.0,26541.9,  8.0,28122.4,  9.0,29780.0,  ... 
10.0,31399.0,  10.4,32092.9,  10.6,32362.7,  11.0,32285.6,  ... 
12.0,32170.0,  15.0,32170.0,  20.0,32170 
AFGEN  NSDATA=  0. 0,2252. 0.  1.0,2245.4,  2.0,2238.8,  3.0,2232.2,  ... 
4.0,2225.6,  5.0,2278.4,  6.0,2436.8,  7.0,256e.8,  ... 
7.7,2621.6,  8.0,2611.0,  9.0,2579.4,  10.0,2568.8,  ... 
11.0,2542.4,  13.0,2516.0,  15.0,2516.0,  20.0,2516 
AFGEN  MFDATA=  0. 0, 78. 8,  1.0,78.8,  3.0,78.8,  4.0,87.1,  ... 

4.5,95.3,  5.0,103.6,  6.0,111.9,  7.0,128.5,  8.0,136.8,  ... 
9.0,161.7,  10.0,174.1,  10.4,178.3,  11.0,170.0,  ... 

15.0,170.0,  20.0,170.0 


AFGEN  QDDATA=  0. 0, 25. 0,  2.0,25.0,  4.0,25.0,  5.0,25.0,  ... 

6.0,49.6,  7.0,92.6,  8.0,129.4,  9.0,178.6,  ... 

10.0. "33. 9,  11.0,264.6,  12.0,252.3,  13.0,246.1,  ... 
15.0,240.0,  20.0,240.0 


INITIAL 

* 

*  ESTABLISH  INITIAL  CONDITIONS. 

NG I =20220. 0 
NS  I =2252 . 0 
MFI  =  78. 75 
QDI  =  25. 0 


*  SET  INITIAL  STATE  PERTURBATION  TO  ZERO 

* 

DNG  =  0.  0 
DNS  =0.0 
DE  =  0.  0 
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DYNAMIC 


NOB  =  26000.0 

NSB  =  1500.0 

MFB  =  122. 5 

QDB  =  220. 0 

QCB  =  80.  0 

JOB  =  ( NGB/QCB) * JG 

JDB  =  ( NSB/QDB) * JD 


DATA  CURVE  FORMULATION 

NGD  =  AFGEN(  NGDATA, TIME ) 

NSD  =  AFGEN( NSDATA, TIME ) 

MFD  =  AFGEN( MFDATA .TIME) 

QDD  =  AFGEN( QDDATA , TIME ) 

MFDLY  =  TRANSP(  100, MF 1 , 0.  70 , MFD) 

STATE  SPACE  LINEAR  MODEL  FORMULATION 

All  =  (20. 282-32. 560*( NGL/NGB ))/(  JGB ) 

A12  =  ( -0. 448 )/( JGB) 

A13  =  (  13. 624 )/(  JGB) 

A21  =  ( 2.  1131  +  . 91786*( NGL/NGB ) )/(  JDB ) 
A22  =  (-. 0998-0. 445*( NGL/NGB) )/( JDB) 
A23  =  ( 18851 )/( JDB) 

A3 1  =  0. 0 


A32  =  0.  0 


A3 3  =  -1. O/TC 

Bll  =  0. 0 

B12  =0.0 

B21  =. 035250/(  JDB) 

B22  =  -1. 0/( JDB) 

B31  =  1. O/TC 


B32  =  0. 0 


DERIVATIVE 

NOSORT 


COMPUTE  INPUT  TO  THE  NONLINEAR  MODEL,  MF( T ) ,  WW( T ) . 

RUN  #7 

DMF  =( MFDLY  -  MF I ) /MFB 
DQD  =(  QDD  -  QDI)/QDE 

DNGDOT  =  A11*DNG  +  A12*DNS  *  A13*DE  +  E11*DMF 

DNSDOT  =  A21*DNG  +  A22*DNS  *  A23*DE  *  E2 1 *DMF  +  E22*DQD 

DEDOT  =  A33-DE  +  E31*D"F 


DYNAMIC  EQUATIONS  FOP  LINEAR  MODEL. 
DNG=INTGF.L(  0.  0,  DNGDOT) 
DNS=INTGRL( 0.  0, DNSDOT) 

DE  =IHTGRL( 0.  0. DEDOT) 

NGL  =  NGI  *  PNG* NOB 
NSL  =  MSI  *  D1.'S*NSB 
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*  OF  'NG',  AND  'NS*  AS  CALCULATED  BY  THE  NONLINEAR  AND  STATE  SPACE 

*  MODELS.  THE  STATEMENTS  BELOW  COMPUTE  THE  VALUES  OF  ' NG1 ,  AND  'NS' 

*  AS  RECORDED  FROM  GAS  TURBINE  TEST  DATA. 

* 

TERMINAL 
METHOD  RK5 

CONTROL  FINTIM=20.  0,DELT=0.  001 
PRINT  0. 25 , NGD, NGL, NSD , NS L, MFD, ODD 
END 
STOP 
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*  * 

*  * 

*  BOEING  MODEL  502-6A  CAS  TURBINE  * 

*  * 

*  DYNAMIC  COMPUTER  SIMULATION  * 

*  * 

*  MODIFIED  BY  * 

*  S. D.  METZ  * 

*  FROM  ROUTINES  BY  * 

*  V.  J.  HERDA  AND  V.  A.  STAMMETTI  * 

*  * 


*  THIS  PROGRAM  SIMULATES  THE  DYNAMIC  RESPONSE  OF  THE  NPS  * 

*  BOEING  GAS  TURBINE  TEST  FACILITY  USING  A  MULTILPLE  * 

*  LINEARIZATION  TECHNIQUE.  * 

*  * 


C 

C 

PARAM  JG=0. 009525,  JD=0. 3000,  PI=3. 14159,  TC=1.  00 
* 

*  THE  FOLLOWING  VALUES  LISTED  ON  THE  FUNCTION  CARD  ARE  FOR  FUEL  FLOW, 

*  GAS  GENERATOR  SPEED,  TORQUE  AND  DYNO  SPEED  AS  A  FUNCTION  OF  TIME. 

*  THESE  VALUES  WERE  OBTAINED  FROM  STRIP  CHART  RECORDS  AND  ARE  ENTERED 

*  IN  THE  FORM  ( E. G.  FUEL  FLOW)  .  .  .  TIME( SEC ) ,  FUEL  FLOW . 


C 

C  THIS  SET  IS  FOR  EXPERIMENTAL  RUN  #  9. 

C 

C 

AFGEN  NGDATA=  0. 0,32150. 0,  5.0,32150.0,  10.0,32150.0,  15.0,32150.0,  ... 

20.0,32150.0,  25.0,32150.0,  50.0,32150.0 
AFGEN  NSDATA=  0. 0,510. 0,  2.0,510.0,  3.0,576.1,  4.0,642.14,  ... 

5.0,708.2,  6.0,774.3,  7.0,840.4,  8.0,922.9,  9.0,1005.5,  ... 
10.0,1137.7,  11.0,1269.8,  12.0,1435.0,  13.0,1583.6,  ... 

14.0,1798.4.  15.0,1963.5,  16.0,2128.7,  17.0,2227.8,  ... 

18-0,2326.9,  19.0,2369.0,  20.0,2409.5,  2 1.0,2459.1,  ... 

22.0,24s2.1,  23.0,2508.6,  24.0,2525.0,  25.0,2525.0,  ... 

50. 0, 2525. 0 

AFGEN  MFDATA=  0.  0, 155.  8,  3.0,155.8,  5.0,155.8,  7.0,160.0,  ... 

10.0,160.0,  13.0,164.1,  15.0,164.1,  18.0,164.1,  ... 

20.0,168.3,  23.0,168.3,  24.0,168.3,  25.0,168.3,  50.0,168.3 


AFGEN  QDDATA=  0.  0, 394.  0.  2.0,394.0,  3.0,382.5,  4.0,371.0,  ... 
5.0,359.5,  7.0,336.5,  8.0,325.0,  9.0,313.50,  ... 

10.0,302.0,  11.0,284.8,  12.0,261.8,  13.0,250.3,  ... 

14.0,233.0,  15.0,221.5,  16.0,221.5,  19.0,221.5,  ... 

20.0,221.5,  21.0,210.0,  23.0,210.0,  25.0,210.0,  50.0,210.0 


INITIAL 


*  ESTABLISH  INITIAL  CONDITIONS. 

* 

NG I =32 1 50.  0 
NSI  =  510.  0 
MFI  =  155.8 
QDI  =  394.  0 


*  SET  INITIAL  STATE  PERTURBATION  TO  ZERO 


DNG  =  0.0 
DNS  =  0.0 
DE  =0.0 
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NGB  =  26000. 0 

NSB  =  1500. 0 

MFB  =  122. 5 

QDB  =  220. 0 

QCB  =  80. 0 

JGB  =  ( NGB/QCB ) * JG 

JDB  =  ( NSB/QDB ) * JD 

DYNAMIC 

* 

*  DATA  CURVE  FORMULATION 


NGD  =  AFGEN( NGDATA , T I ME ) 

NSD  =  AFGEN(NSDATA,TIME) 

MFD  =  AFGEN( MFDATA, TIME ) 

QDD  =  AFGEN( QDDATA, TIME ) 

MFDLY  =  TRANSP(  100, MFI , 0.  70, MFD) 

STATE  SPACE  LINEAR  MODEL  FORMULATION 


All 

— 

(20. 282-32. 560*(NGL/NGB) )/(  1GB) 

A12 

= 

( -0. 4480 )/(  JGB) 

A13 

= 

(  13.  624|l/<  JGB) 

A21 

= 

(2. 1131*. 91786*( NGL/NGB) )/( JDB) 

A22 

= 

( -. 0998-0. 445*( NGL/NGB) ) /(  JDB) 

A23 

=  ( 

-0. 18851 )/( JDB) 

A3 1 

= 

0.  0 

A32 

= 

0.  0 

A33 

= 

-1. O/TC 

Bll 

= 

0.  0 

B12 

= 

0.  0 

B21 

=  . 

035250/( JDB) 

B22 

= 

-1. 0/( JDB) 

B31 

= 

1 .  O/TC 

B32 

= 

0.  0 

DERIVATIVE 

NOSORT 

* 

*  COMPUTE  INPUT  TO  THE  NONLINEAR  MODEL,  MF{ T ) ,  WW(T). 

*  RUN  #9 

DMF  =( MFDLY  -  ME I )/MFB 
DQD  =( QDD  -  QDI )/QDB 

DNCDOT  =  All* DUG  ♦  A12*DNS  +  A13*DE  +  E11*DMF 

DNSDOT  =  A2 1 *DNG  +  A22*DNS  +  A23*DE  +  B7 1  *  DMF  +  B22*DQD 

DEDOT  =  A33*DE  +  B3 1 *DMF 


f 


*  DYNAMIC  EQUATIONS  FOP  LINEAR  MODEL. 
DNG=INTCRL( 0.  0 , DNCDOT ) 
DNS=INTGFL(  0  O.DNSDOT) 

DE  =INTGRL( O. C, DEDOT) 

NGL  =  NG I  *  DUG* NGB 
NSL  =  NS  I  +  DNS  *NSB 

SORT 
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* 

*  THE  STATEMENTS  IN  THE  PREVIOUS  (DERIVATIVE)  SECTION  YIELD  VALUES 

*  OF  ' KG',  AND  'NS'  AS  CALCULATED  BY  THE  NONLINEAR  AND  STATE  SPACE 

*  MODELS.  THE  STATEMENTS  BELOW  COMPUTE  THE  VALUES  OF  'NG' ,  AND  'NS' 

*  AS  RECORDED  FROM  GAS  TURBINE  TEST  DATA. 

* 

TERMINAL 
METHOD  RK5 

CONTROL  FINTIM=50.  0,DELT=0.  001 

v  PRINT  0. 5 , NGD,NGL, NSD, NSL, MFD, QDD, DNSDOT, A21 , DNG,  ... 

A22 , DNS , A23,DE,B21, DMF, B22 , DQD 

END 

STOP 

> 
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APPENDIX  E 

SOURCE  DATA  RUNS: 


Run  1:  Data  run  for  dynamic  Ns  run  at  constant  Ng. 

Ng  -  20230  -  20240  rpm  Ns  -  505  -  2215  rpm 

Qd  -  125  -  25  lbft  Mf  -  78.4  -  82.1  Ibm/hr 

Run  3:  Data  run  for  dynamic  Ng  run  at  constant  Ns. 

Ng  —  20220  -  32180  rpm  Ns  --  506  -  513  rpm 

Qd  -  125  -  430  lbft  Mf  -  77.5  -  160.3  lbm/hr 

Run  5:  Data  run  for  dynamic  Ng  run  at  constant  Ns 

Ng  -  20220  -  32130  rpm  Ns  -  1492  -  1510  rpm 

Qd  -  71  -  144  lbft  Mf  -  83.8  -  165  lbm/hr 

Run  7  :  Data  run  for  dynamic  Ng  run  at  Ns  (2252-2516  rpm). 

Ng  -  20220  -  32170  rpm  Ns  -  2252  -  2516  rpm 

Qd  -  25  -  240  lbft  Mf  -  78.8  -  170  lbm/hr 

Run  9:  Data  run  for  dynamic  Ns  run  at  constant  Ng. 

Ng-  32150  rpm  Ns  -  510  -  2525  rpm 

Qd  -  394  -  210  lbft  Mf  -  155.8  -  168.3  lbm/hr 
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